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ABSTRACT 

This  report  addresses  the  problem  of  estimating  the  parameters  of 
superimposed  signals  observed  by  an  array  of  sensors.  Some  of  the  proposed 
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noise.  The  methods  used  are  direct  iterative  maximum  likelihood,  the  EM 
algorithm,  the  e igenstructure  approach  and  the  polynomial  approach.  In 
addition  to  the  traditional  estimation  of  the  source  locations  we  also 
address  the  estimation  of  parameters  related  to  the  radiation  patterns  of 
the  sources. 
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I.  INrR(M)UCT10N 

This  report  addresses  the  problem  of  estimating  the  parameters  of 
super iiiQ>osed  signals,  occurring  in  a  variety  of  fields  ranging  from  radar, 
sonar  and  oceanography  to  seismology  and  radio-astronomy. 

In  this  section  we  formulate  the  superimposed  signals  problem,  review 
the  relevant  literature,  and  present  a  summary  of  the  content  and  the 
contributions  of  this  report, 

1.1  Formulation  of  the  problem 

Our  formulation  of  the  superin^osed  signals  problem  is  motivated  by  the 
specific  problem  of  localizing  N  radiating  sources  using  an  array  of  M 
sensors.  The  signal  at  the  output  of  the  m-th  sensor  can  be  described  by 


X 

m 


(t) 


N 


O  S  ( t-T  ) 

mn  n  mn 


n=l 


V  (t)j  m=l,2,..,,M. 


-T/2  <  t  <  T/2 


(1) 


where  {Sj^(t)}^_]^  are  the  radiated  signals,  additive  noise 

processes,  and  T  is  the  observation  interval.  The  intensities  and  the 

delays  are  parameters  related  to  the  directional  patterns  and  relative 

locations  of  the  n-th  source  and  the  m-th  sensor.  Note  that  a__  is  a 
function  of  the  radiation  pattern  of  the  source  in  the  direction  of  the 
sensor,  of  the  radition  pattern  of  the  sensor  in  the  direction  of  the  source 
and  of  the  distance  between  the  source  and  the  sensor.  However,  is  only 

a  function  of  the  distance  between  the  source  and  the  sensor.  Hence  the 
estimation  of  {a„„}  and  yields  important  information  on  the  locations 


3 


and  radiation  patterns  of  tlie  sources. 

A  convenient  separation  of  the  parameters  of  interest  is  obtained  by 
using  Fourier  coefficients  defined  by 


X  (» 
m 


*  iT*  V  nr  /<• 


X  (t)  e 


-joj  .t 


T  ''-T/2 


m 


dt. 


where  =  2n(i]^+i)/T,  i  =  l,2,...,Ij  and  ij^  is  a  constant.  In  principle 
the  number  of  required  coefficients  tends  to  infinity.  However,  since  we 
consider  only  finite  bandwidth  signals,  we  can  use  only  !<“'  coefficients. 
Taking  the  Fourier  coefficients  of  (1)  we  obtain: 


X  (w.) 
m  1 


N 


-JW  .T 

a  e  ^  ”*^8  (w.)  +  V  (w.), 
mn  nr  mi 


n=l 


(2) 


where  Sj^(<))j^)  and  are  the  Fourier  coefficients  of  Sjj(t)  and  Vjj,(t) 

respectively.  Equation  (2)  may  be  expressed  using  vector  notation  as 
follows: 


X(a)^)  =  A(w.)S(w^)  +  V(a)^)j 


i  =  1,2, ...,I; 


(3) 


where 
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X(a)^)  =  [X^(uJ,  X2(a)^),...,Xjj((o.)]^ 
S(a)^)  =  [Sj(<D^),  S2((i>j),...,Sj^(a>j)]^ 

V(a),)  =  [V-(a),),  V,(a). V„(<oJ]^ 

~  1  lx  21  Ml 

A(<1).)  =  [a  (e  ),  a  a. (§„)] 

X  1  **  1  1  "“1  “"XN 


Si<«n> 


“2® 


*  •  •  • » ®  J 

ffln 


We  use  to  represent  all  the  parameters  of  interest  associated  with  the  n- 
th  signal,  namely  and  .  Our  main  goal  is  to  estimate  the 

set  t§n}n=i»  Note  that  if  the  spectrum  of  the  signals  is  concentrated 

around  (Oj^,  with  a  bandwidth  that  is  small  compared  to  2it/T,  then  (3)  reduces 
to  a  single  relation  between  the  observation  vector  and  the 

parameters,  i.e.  1=1,  In  this  case  it  is  customary  to  use  many  short 

observation  intervals  or  siinply  time  samples,  and  the  model  becomes: 

X(j)  =  As(j)  +  V(j);  j=l,2,...,Jj  (4) 

where  the  dependence  on  the  single  frequency  is  suppressed,  and  j  is  the 
index  of  the  different  sauries.  Note  that  the  main  difference  between  the 
narrowband  case  and  the  wideband  case  is  that  A  is  the  same  in  all  the  J 

equations  specified  by  (4)  while  ACw^)  is  different  in  each  of  the  I 

equations  given  by  (3).  In  this  report  we  concentrate  on  the  wideband  case 


whenever  the  proposed  procedure  can  handle  both  the  wideband  case  and  the 
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narrowband  case. 

Under  tbe  assuo^tion  tbat  tbe  number  of  sources  is  known,  the  least 
squares  estimates  of  {0^^}  is  given  by: 


{0  }  =  arg  min  Qj 

“  {e  }ee 

-n 


I 

Q=  5  ||X(Wi)  -  A((0i)S(a)i)  11^ 
i=l 


(5) 


where  ll'll  denotes  the  Euclidean  norm  and  ©  is  the  given  parameter  space. 

Equation  (5)  also  represents  the  maximum  likelihood  estimates  under  the 

assumption  that  the  noise  vectors  {VCw^)}  are  i.i.d.  zero  mean  Gaussian  with 
2 

covariance  cr  I. 

The  minimization  required  in  (5)  is  not  trivial  since  the  vectors  SCw^) 
and  the  matrix  A((i)£)  are  not  known  to  the  observer.  However,  whenever 
is  known  Q  is  minimized  by  choosing 

|(w^)  =  [A(u).)®A(u).)]"^A(a).)^(yp  (g) 


as  the  estimates  of  S(^,£)^  where  (  )®  denotes  the  Hermit ian-transpose 

operation.  Substituting  (6)  in  (5)  we  obtain 

I 

{§  }  =  arg  max  3  X(a)  )®A<<o  J  [A(<d.)\(w.)]"^A((o  .)\(w  •)  .  (7) 

{©n]8©  i=i 


the  maximization  in  (7)  requires  a  multidimensional  search  over  all  the 
parameters  and  and  since  this  problem  is  difficult  many  papers  and 
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books  bave  proposed  subopt imam  estimation  schemes. 

1.2  Literature  Survey 

A  comprehensive  literature  survey,  including  more  than  120  references 
is  included  in  [1].  Also  see  [2]  for  many  other  references  not  discussed  in 
[1] .  It  is  beyond  the  scope  of  this  report  to  describe  every  algorithm  and 
every  set  of  assun^tions  in  the  hundreds  of  estimation  schemes  that  have 
been  proposed  until  now.  Instead  we  confine  our  attention  to  the  techniques 
that  are  currently  the  most  promising. 

The  EM  Algorithm  [6]-[8] 

The  EM  (Expectation-Maximization)  algorithm  was  recently  proposed  by  M. 
Feder  and  E.  Weinstein  for  treating  the  problem  described  above.  The  EM 
procedure  is  essentially  an  iterative  algorithm  that  is  guaranteed  to 
converge  to  a  stationary  point  of  the  likelihood  function.  Hence,  if  it 
converges  to  the  global  maximum  of  the  likelihood  function  the  estimates  are 
exactly  the  maximum  likelihood  estimates.  The  main  disadvantages  of  this 
algorithm  are: 

1)  If  the  likelihood  function  is  not  unimodal  the  algorithm  may 
converge  to  a  local  maximum.  Hence,  it  may  be  necessary  to 
overcome  this  problem  by  appropriate  measures. 

2)  The  algorithm  is  often  slow  to  converge  and  the  amount  of 
computation  required  for  each  iteration  may  be  large. 

In  [6]-[8]  Feder  and  Weinstein  derived  the  EM  algorithm  for  the  general 
linear  Gaussian  case  for  known  signals  in  noise  and  random  signals  in  noise 
with  known  Gaussian  statistics.  Here  we  extend  these  results  to  the  more 
realistic  case  of  non-random  unknown  signals  in  noise.  We  also  make  use  of 
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the  EM  algorithm  in  order  to  apply  the  polynomial  approach  to  nonuniform 
arrays . 

The  Covariance  Eigenstrncture  Approach 

The  covariance  eigenstructure  approach  was  first  proposed  in  the  time- 
domain  by  Schmidt  [3],  who  called  it  the  Multiple  Signal  Identification  and 
Classification  (MUSIC)  method,  and  by  Bienvenu  [10]  who  developed  an 
equivalent  frequencjr-domain  procedure.  Since  its  introduction,  a  large 
number  of  extensions  and  refinements  of  this  method  have  been  proposed,  and 
this  technique  is  therefore  considered  as  one  of  the  most  practical  for 
solving  the  super inq>osed  signal  identification  and  retrieval  problem.  Note 
however  that  since  this  method  does  not  atten^t  directly  or  indirectly  to 
maximize  the  likelihood  function,  it  is  suboptimal.  Yet,  it  yields  good 
results  for  sufficiently  high  signal  to  noise  ratio  (SNR).  Its  main 
advantages  are: 

1)  It  relies  on  an  algorithm  which  is  not  iterative,  and  hence  it 
eliminates  the  problem  of  converging  to  local  stationary  points. 

2)  The  amount  of  coiiq>utation  is  less  than  the  amount  required  for  the 
EM  algorithm. 

Its  main  drawbacks  are: 

1)  The  algorithm  cannot  be  used  in  problems  in  which  there  are  only  a 
small  number  of  observations. 

2)  There  is  no  natural  way  to  extend  the  algorithm  to  handle  wideband 
signals.  For  different  extensions  for  wideband  signals,  see  [91- 
[15]. 

In  this  report  we  show  how  the  MUSIC  algorithm  can  be  used  in  the 
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rather  interesting  and  practical  case  of  non-omnidirectional  sources  in  the 
near  field  of  the  array. 

The  Polynomial  Approach  [5],  [16]-[18]« 

The  polynomial  approach  in  the  context  of  maximum  likelihood  (BfL) 
estimation  was  introduced  only  recently  by  Bresler  and  Macovski  [5].  This 
approach  is  limited  to  the  special,  but  in^ortant,  case  of  linear,  uniformly 
spaced,  narrowband  arrays. 

Although  the  starting  point  of  this  technique  is  precisely  the  ML 
estimation  problem  described  above,  the  algorithm  proposed  in  [5]  is  not 
guaranteed  to  yield  results  that  are  confined  to  the  a-priori  known 
parameter  space.  However,  in  simulations,  the  algorithm  converges  within  5- 
10  iterations  to  the  right  results,  for  high  enough  SNR.  In  this  report  we 
extend  the  polynomial  approach  to  nonuniform  arrays. 

Summary  of  Content  and  Contributions 

In  Section  II  we  first  briefly  review  the  EM  algorithm  and  then 

following  the  approach  proposed  by  Feder  and  Weinstein,  we  derive  an  EM 
algorithm  for  the  general  case  of  superiiiq)Osed,  unknown  deterministic 
signals  in  noise.  This  extension  of  the  work  in  [6]-[8]  is  important  for 

obvious  reasons.  However,  the  results  are  very  general  and  require  some 

refinement  for  practical  coiiq)utation.  This  is  done  in  the  last  part  of 
Section  II  in  which  we  show  how  to  obtain  an  estimation  of  the  location  and 
intensity  parameters  by  a  relatively  efficient  EM  procedure.  In  Section  III 
we  describe  a  novel  and  efficient  algorithm  for  coitq>uting  the  ML  estimates 
of  superiit^osed  signals.  The  algorithm  is  equally  applicable  to  wideband 
sources  and  narrowband  sources  and  does  not  require  a  knowledge  of  the 
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statistical  properties  of  the  signals.  Typically,  it  requires  less 
iterations  than  the  EM  algorithm.  In  Section  IV  we  briefly  review  the 
eigenstructure  approach  and  then  we  employ  some  of  the  ideas  of  Section  II 
to  extend  this  approach  to  the  estimation  of  the  radiation  patterns  of  the 
sources  as  well  as  the  location  of  the  sources.  Some  users  may  prefer  this 
suboptimal  approach,  since  it  is  very  fast  and  yields  good  results  at  high 
enough  SNR.  Section  V  is  devoted  to  linear,  narrowband  arrays  with 
nonuniform  sensor  spacings.  We  make  use  of  both  the  polynomial  approach  and 
the  EM  algorithm. 
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II.  APPLICATION  OF  THE  EM  ALGORITHM  TO  THE  ESTIMATION  OF  SUPERIMPOSED 
UNKNOWN  SIGNALS  IN  NOISE 

In  this  section  we  briefly  review  the  EM  method  and  then  apply  it  to 
the  ML  problem  described  above. 

II. 1  The  EM  Method 

Let  X  denote  the  observation  vector  and  0  represent  the  parameter 
vector.  If  fj.(X|0)  is  the  conditional  probability  density  function  of  x 
given  0,  then  the  ML  estimate  of  6  is: 

0  =  arg  max  f  (Xj0)  =  arg  max  ln{f  (X|0)},  (8) 

0e©  -  080  - 

where  0  is  the  parameter  space. 

In  many  cases  of  interest  one  would  like  to  observe  Y,  the  •coiiq)lete 
data",  instead  of  X,  the  "incon^lete  data",  where  the  relation  between  X  and 
Y  is  given  by  some  non-invert ible  mapping: 

H(Y)  =  X. 

From  Bayes'  rule  we  have 

ln{f^(X|0)}  =  lm{f^(Y|0)}  -  lu{f^j^(Y|X,0)}  (9) 

Taking  the  expectation  of  (9)  over  y  given  X  and  under  the  assumption  that 
the  parameter  is  equal  to  0',  we  obtain 
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L(e)  =  in{f^(x|e)i  =  Q(e|e')  -  H(e|e'), 


where 


Q(e|e')  =  E{in{f^(Y|e)}|x,e'} 
H(eje')  =  E{in{f  I  (i|x,e)} |x,e'}. 


Using  Jensen's  inequality  it  is  easy  to  verify  that 

H(eje')  <  H(e' j©'). 

The  EM  procedure  may  be  described  by  the  following  sequence  [20] : 

(a)  Initialization:  set  p=0,  and  =  ©q, 

(b)  E-step:  Determine  Q(©j©^P^). 

(c)  M-step:  Choose  ©^P+D  to  be  the  value  of  ©e©  that 

Q(©|©<P)). 

(d)  Check  the  convergence  of  ©.  No  -  p=p+l,  go  to  (b). 

Yes  -  stop. 

In  every  cycle  of  the  algorithm  the  likelihood  function  L(©)  is 
s inc  e : 

L(0(P+1))  1©^^^)  -  1©^^^) 


(10) 


(11) 


maximizes 


inc reased. 
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The  inequality  holds  due  to  (11)  and  due  to  the  M-step.  For  proof  of 
convergence  of  this  procedure,  see  [21]. 


II. 2  Application  to  Array  Processing 

In  this  section  we  apply  the  EM  method  to  the  case  of  super in^osed 
unknown  nonrandom  signals.  We  concentrate  on  the  wideband  case  described  by 
equations  (5).  The  modification  for  the  narrowband  case  is  straightforward. 

Following  [6]-[8]  we  choose  the  "complete  data"  as  the  observation  of 
each  of  the  signals  separately.  Hence 


a.(e  )S  (w.)  +  V  (m.), 
-1  -n  n  1  -n  i 


(12) 


where  the  fictitious  noises  chosen  to  be  mutually  uncorrelated, 

zero  mean  Gaussian  vectors  with  covariance  o^I  satisfying 

N 

NorJ  =  and  ^  =  YCw^),  i  =  1,2....,I.  (13) 

n=l 


The  coiiq)lete  data  vector  is 

Y  =  [y'’'^((o^),  Y^(w2),...,y'’^(o>j)]^, 


where 
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Y(a).)  =  [YjU.),  Y^Co). yJ<<o.)]^. 

The  "incoE^lete  data*  (i.e.  the  observed  data)  X  is  defined  by 
X  =  [X^(a)j),  X^Uj) . 


and  is  obtained  from  Y  by  the  linear  transformation 


5  =  GX. 

where  G  is  a  block  diagonal  matrix  with  I  blocks: 

■e  0  ■ 

r  -  ® 

G  —  •  > 

.0  ’H  . 

and  the  matrix  H  is  constructed  of  N  identity  matrices  of  size  MxM, 

“  -  'Vm  •••  V  • 

We  now  turn  to  evaluate  the  functional  G(GjO').  We  recall  that  Y  is 
Gaussian  with  known  covariance  a|l  and  unknown  mean  ii(§)  (where  the 
parameter  vector  G  includes  not  only  the  parameters  {0^,,,}  and  but  also 

the  signal  parameters  Sj^Cto^)  with  llill,  llnlN) ,  hence 


14 


ln{f  (Y|X,e)}  =  -In  det(nGjl)  - -yj  |X  - 

“i 

and 


Q(e|e»)  =  K  -  -4-  Ilf  -  ii(S)|P 

•^1 

where  E  represents  terms  independent  of  6  and 

t  =  Etll?,  0'}  =  li(e')  +  G°(GG®)"^(X  -  Gii(e')). 

Equation  (15)  may  be  rewritten  as 

N  I 

9<«lr)  -  K  -  4-  I  1 

®1  n=l  i=l 

and  using  the  block  diagonal  structure  of  G,  equation  (16)  becomes: 

Y  (o)  )  =  a  (G'  )Sj;(<o,)  +4  tX(a),)  -  A' (m .) S' (w . ) ] 

— n  1  — 1  —  nni  N—  i  i—  i 

The  proposed  EM  algorithm  may  be  summarized  as  follows: 

(a)  Guess  initial  values  for  the  parameters  fa,  and  {-c 
construct  the  matrices  A(<Oj),  i=l,2,...,I.  Confute 

estimates  for  S((i)£)  using  (6). 


(14) 


(15) 


(16) 


(17) 


(18) 


mnJ ' 
initial 


(b)  E-step:  Substitute  in  (18)  the  current  estimates  of  the  parameters 
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and  compute  for  llnlN,  llill. 

(c)  M-step;  Find  the  maximum  of  (17)  for  each  0^^.  This  is  simply: 


e(p+i)  ^ 
-n 


arg  max 
0 

-n 


g(p+l) 

n 


(m.) 

1 


aH(0(p+l)  • 

“1  —n  i  •  •—1  n 


)ir 


(d)  Check  convergence  of  {0jj} .  If  not:  go  to  step  (b). 

If  yes:  done. 

Observe  that  the  EM  algorithm  presented  here  solves  iteratively  the 
original  maximization  problem  over  the  parameters  of  N  signals.  At  each 
iteration  we  have  to  solve  N  reduced  maximization  problems,  one  for  each 
signal.  However,  even  the  reduced  maximization  problem  specified  by  the  M- 
step  is  not  trivial,  since  2(M-1)  parameters  {g^^} ,  {Tjpjj}  are  not  known. 
Therefore,  in  the  following  sections  we  will  attempt  to  reduce  further  the 
computational  reguirements  of  each  iteration. 


II. 3  Further  Sin^lif ication  of  the  EM  Procedure 

In  this  section  we  show  how  the  parameters  of  each  source  can  be 
estimated  with  minimal  effort.  We  first  define  two  vectors: 


a 

-n 


‘  <“l.-  “2»-- 


.,a  )■ 

mn 


in  ■  <’^ln-  ■'2n' 


.T  )■ 
mn 


which  we  call  the  intensity  vector  and  the  delay  vector  of  the  n-th  source. 
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Note  that  one  can  always  choose  the  first  con^onent  of  to  be  0,  and 
Ibnil  =  1  without  loss  of  generality.  This  is  true  since  we  have  extra 
degrees  of  freedom  due  to  the  estimation  of  both  and  The  n-th 
column  of  ACw^)  is  a  function  of  both  and  however,  these  parameters 
separate  as  follows: 


(19) 


where  is  a  diagonal  matrix  defined  by 


T  ix  )  =  diagCl,  e 

1  — -n 


mn. 

» •  •  •  f  0  )  ^ 


Now,  the  maximization  problem  in  the  M-step  becomes: 


(p+1)  _  (p+1)  _(p+l) 

=  o  ,  X 

-n  -n  -n 


arg  max  }  a®(e^)YfP^(o  J  )®a  .  (6  ) 

-n  i-l”  ~  ^ 


=  arg  max  At  )YfP\a)  J  ) )®r, (t  )]a  . 

_  — n  4.  1  — n  — n  i  — n  i  i  — n  -n 


t  ,o  .  . 

-n  -n  1=1 


The  solution  of  the  above  maximization  problem  is  given  by: 
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,(p+l) 


n 


=  arg  max 
X 

-n 


-  ID&X  r  «  r  V  ^ 

X  {C(t  )} 
"“H 


^(p+l)  _  ^max 
“h  — 


(20a) 


(20b) 


vrhere  X™®^{C(Tjj^) }  is  the  largest  eigenvalue  of  the  matrix  Cix^^)  defined  by 
I 

C(t  )  =  Re{  )  At  )y^P^a),)(YyNa).))V(t  )} 

— n  i.  t  — n  — n  i  — n  i  i  — n 

i=l 

and  is  the  associated  eigenvector. 

Equation  (20)  requires  a  search  over  (M-1)  parameters  (the  co]:q>onents 
•  However,  even  this  problem  can  be  reduced  by  recalling  that  the 
delay  parameters  are  not  independent.  One  can  express  each  of  the  delays  in 
tjj  as  a  function  of  only  two  or  three  source  location  parameters.  This  is 
true  since  the  delays  are  only  a  function  of  the  distance  between  the  source 
and  the  sensor  (we  assume  that  the  speed  of  propagation  in  the  medium  is 
known  and  that  the  sensor  location  is  known).  Now,  the  search  is  limited  to 
a  three  dimensional  search  over  all  possible  individual  source  locations. 
If  one  is  interested  in  the  planar  case  or  azimuth  only  system,  the  search 
is  confined  to  only  two  or  one  dimension,  respectively. 

Finally,  we  note  that  the  above  method  provides  a  very  useful  tool  for 
estimating  the  vectors  which  in  turn  provide  valuable  information 

regarding  the  directional  properties  of  the  sources  and/or  the  sensors,  and 
also  might  be  used  to  evaluate  the  attenuation  of  the  medium  in  various 
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directions. 
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III.  Direct  Maximam  Likelihood  Approach 

In  this  section  we  present  a  novel  and  efficient  algorithm  for 
computing  the  maximum  likelihood  estimates  of  multiple  signals  observed  by 
an  array  of  sensors.  The  algorithm  provides  estimates  of  parameters  related 
to  the  directional  patterns  of  the  sources  {0^^}  as  well  as  estimates  of  the 
location  parameters  of  the  sources  .  Furthermore,  the  algorithm  is 

equally  applicable  to  wideband  sources  and  narrowband  sources  and  does  not 
require  a  knowledge  of  the  statistical  properties  of  the  signals.  In  this 
section  we  concentrate  on  the  wideband  case.  The  modification  for  the 
narrowband  case  is  straightforward,  and  can  be  found  in  LIDS-P-1670. 

We  basically  want  to  find  a  solution  for  equation  (5).  Relation  (6) 

* 

enables  us  to  update  the  estimates  S(<i>£)  whenever  we  have  new  estimates  for 
A(<0j^).  The  main  principle  of  the  algorithm  is  to  perform  successive 
minimization  operations  on  the  parameters  of  each  signal,  holding  all  the 
rest  of  the  parameters  fixed.  For  example,  suppose  that  we  want  to  perform 
a  minimization  with  respect  to  the  k-th  signal  parameters,  then  Q  can  be 
rewritten  as 

I 

Q=  5  11^(0).)  -  a^(0j^)Sj^(a,.)||2  (21) 

i=l 

where  aj(6^)  is  the  k-th  column  of  ACw^),  is  the  k-th  component  of 

and  is  given  by 
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Y^(w.)  =  X(w^)  -  A(<o  (22) 

where  S^((Oj^)  is  sinqdy  SCw^)  with  the  k-th  con^onent  replaced  by  zero. 

The  minimization  of  (21)  with  respect  to  Ok,  using  (6)  with  A(w£) 
replaced  by  aj(0j,),  is  given  by 

I 

0^  =  arg  min  ^  |  .)  -  a.  (Oj^)  [a®(ej^)a.  (ej^)]"V(ej.)Y^(u)  )  1 

0,  80  .  , 

-k  1=1 

I 

=  arg  mar  J  I  1^/ I  UjCOjj)  I  (23) 

5l'®  i-l 

We  now  apply  the  assun^tion  |  |a j|^(0jj.)  ( P  =  1  and  the  decomposition  (19)  to 
obtain: 

=  arg  max  (24a) 

a  = 

— k  — 

where  is  the  largest  eigenvalue  of  the  matrix  given  by: 

I 

R^  =  Re{  ^  (25) 

i=l 

and  is  the  associated  normalized  eigenvector. 
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The  maximization  described  by  (24a)  can  be  performed  by  a  siii^>le  search 
over  the  space  of  induced  by  all  possible  individual  source  locations, 
or  by  a  single  gradient  subalgorithm. 

The  algorithm  is  summarized  as  follows: 


(a) 

Initialization:  Select 

,  Set  k=l. 

(b) 

Confute  S(w£)  according  to  (6). 

(c) 

Confute  Kj.  according  to  (22)  and 

(25). 

(d) 

Find  ^  according  to  (24). 

(e) 

Update  the  k-th  column  of  A((i)£)  for  llill  with  the  new  set 

k=k+l,  if  k>N  then  k=l. 

(f) 

Check  the  convergence  of 

,  If  yes:  done?  if  no: 

go  to  (b). 

Observe  that  at  each  updating  step  (i.e.,  steps  (b)  and  (e)),  we 
decrease  the  cost  function  Q  defined  in  (5).  Since  0^0  the  algorithm  will 
converge  at  least  to  a  local  minimum  of  Q.  Depending  on  the  initial 
estimates  of  a^,  and  on  the  structure  of  0,  the  local  minimum  may  or  may 

not  coincide  with  the  global  minimum. 

This  algorithm  may  be  viewed  as  a  modification  of  a  special  case  of  the 
EM  algorithm.  According  to  the  theory  of  the  EM  method,  the  estimates 
generated  in  the  M-step  should  be  used  in  the  E-step.  This  may  be  applied 
to  the  present  algorithm  as  follows.  Instead  of  updating  using  (6)  in 

step  (b),  S(<i)j^)  is  updated  by  replacing  only  the  k-th  component  by  the 

estimates,  »  which  can  be  computed  in  step  (d),  following  the 

computation  of  and  8j^.  Note  that  a?(ejj^)Y^(w is  siB5)ly  the  value  of 
Sj.((i)£)  that  minimizes  (21)  whenever  a^(6j^)  is  known.  It  is  clear  that  the 
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last  procedure  typically  will  require  more  iterations  than  the  proposed 
procedure  since  the  updating  of  S((i)^)  is  done  without  using  all  the 
currently  available  information. 
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IV.  EIGENSTRUCTURE  APPROACH  FOR  ARRAY  PROCESSING  WITH  UNKNOWN  INTENSITY 

COEFFICIENTS 

The  eigenstracture  approach  for  array  processing  is  examined  for  the 
general  case  in  which  it  is  required  to  estimate  parameters  related  to  the 
directional  patterns  of  the  sources  as  well  as  parameters  related  to 

the  location  of  the  sources  .  In  recent  years  there  has  been  a  growing 
interest  in  eigenstructure  based  methods,  perhaps  due  to  their  applicability 
to  general  array  configurations  and  due  to  their  sia^licity  and  relative 
efficiency.  A  conqprehensiye  discussion  of  the  method  may  be  found  in  [3], 
while  [1]  contains  a  literature  survey  of  most  of  the  recently  published 
results. 

An  assuiiq>tion  common  to  all  previously  published  contributions  in  this 
area  is  that  any  given  source  is  observed  by  all  the  sensors  with  the  same 
intensity.  This  assumption  is  reasonable  only  if  the  sources  are  in  the 
far-field  of  the  array  and  the  sensors  have  identical  radiation  patterns. 
In  this  report  we  remove  this  rather  restrictive  assuitq>tion  and  thus  extend 
the  applicability  of  the  eigenstructure  approach  to  the  case  of  near-field 
sources  and/or  sensors  with  unknown  radiation  patterns. 

Since  there  are  more  than  one  extension  of  the  eigenstructure  approach 
to  wideband  signals  we  concentrate  here  on  the  narrowband  case.  The 
modification  for  each  of  the  wideband  extensions  described  in  [9]-[15]  is 
straightforward. 

The  following  assun^tions  are  made: 

(a)  The  signals  and  noises  are  stationary  over  the  observation 


interval. 
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(b)  The  number  of  sources  in  known  and  it  is  less  than  the  number  of 
sensors . 

(c)  The  columns  of  A,  in  equation  (4),  are  linearly  independent. 

(d)  The  signals  are  not  conq>letely  correlated. 

(e)  The  noise  covariance  matrix  is  known  except  for  a  multiplicative 
constant  o^. 

Recalling  that  equation  (4)  describes  the  narrowband  case,  the 
correlation  matrices  of  the  signal,  noise  and  observation  vectors  are  given 
by 


Rs  =  E{S  S°} 

\  =  Et?  X®}  =  AR^A®  +  (26) 

TJ 

where  (  )  represents  the  Hermit ian  transpose  operation.  The  following 

theorem  form  the  basis  for  the  eigenstructure  approach. 

Theorem;  Let  and  u^,  k=l,2,...,M  be  the  eigenvalues  and 

corresponding  eigenvectors  of  the  matrix  pencil  (R^,  ^),  with  in 

dec ending  order.  Then, 

%+l  “  %+2  “  •••  ^  • 

2)  Each  of  the  columns  of  A  is  orthogonal  to  the  matrix  U  = 

HN+2’*“*Hm^* 

Proof;  See  [14] . 


This  theorem  suggests  that  reasonable  estimates  of  the  parameters 
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^-n^n=l  obtained  by  first  generating  an  estimate  U  of  U  and  then 
searching  over  all  possible  values  of  0^^  for  vectors  a(6j^)  that  are  nearly 
orthogonal  to  U.  This  may  be  written  as 


0 

-n 


arg  mini |U^ a(0  ) 1 1^ 
Q  II  --nil 

-n 


(27) 


where  ll'll  denotes  the  Euclidean  norm.  Since  there  is  an  extra  degree  of 
freedom,  there  is  no  loss  of  generality  in  assuming  that  |  ja(0jj)  1 1  =  1. 
This  also  eliminates  the  trivial  solution  of  (27).  Note  that  (27)  requires 
a  multidimensional  search  over  the  parameters  and  .  To  overcome 

this  difficulty  we  decon^ose  a(0^)  as  follows: 


a(0  )  =  P(t  ) 
—  — n  ~n 


a 

-n 


where 


-n  “  ^“in*  °2n****'“mn^ 


w/  ^  "j“l^2n  "j“l^mn, 

P(lj^)  =  diag(e  ,  e  »...*«  ) 


-n  "  ^■'in'  '^2n’**‘'‘^mn^  * 


Using  this  notation,  (27)  becomes 
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m 

0 

-n 


=  arg  min 

o  ,x 

-n  -n 


)UU®P(i:  )a  , 
-n  -n  -n  -n 


(28) 


and  hence 


-n 


=  arg  min 
X 

-n 


min,  , 

6  (C(t^)}, 

~n 


a 

-n 


(29a) 


(29b) 


where  6®^“{C(Tj^)}  is  the  smallest  eigenvalue  of  the  matrix  C(Tj^)  given  by 

Cix^)  =  Re{r®(T^)UU®r(T^)},  (30) 

and  is  the  associated  normalized  eigenvector.  Equation  (29)  requires  a 

simple  search  over  the  space  of  vectors  induced  by  all  possible 

individual  source  locations. 

The  proposed  algorithm  may  be  summarized  as  follows: 

(a)  Estimate  the  observation  covariance  matrix: 

J 

•IT  H 

Rx=t  2?u)x(j)''. 

j=i 

(b)  Find  the  M-N  eigenvectors,  corresponding  to  the  smallest  M-N 

eigenvalues  of  the  pencil  (R^,1q),  and  construct  the  matrix: 
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^  ^-N+1’  2n+2' *  *  * ’2m^ • 

(c)  Evaluate,  for  all  possible  source  locations,  tbe  "spatial 
spectrum"  given  by: 


*  \  u  ; 

6““{C(t)} 

where  C(t)  is  defined  by  (30). 

(d)  Select  the  N  highest  peaks  of  P(t).  The  corresponding  values  of  t 
describe  the  source  locations,  and  the  corresponding  eigenvectors 
describe  the  intensity  vectors  . 

Exan^les 

To  illustrate  the  behavior  of  the  algorithm,  let  us  consider  two 
ezac^les: 

Example  1.  Consider  a  uniform  linear  array  of  five  sensors  separated  by 
half  a  wavelength  of  the  actual  narrowband  source  signals.  The  sources  are 
two  narrowband  emitters  located  in  the  farfield  of  the  array.  In  this  case, 
if  Yjj  denotes  the  bearing  of  the  n-th  source,  n=l,2,  relative  to  the 
perpendicular  to  the  array  baseline,  the  differential  delay  is  given  by 
=  (m-l)nsin(yjj^) .  The  first  source  at  a  bearing  of  -9  degrees  was  observed 
with  the  intensity  vector  aj  =  [1,1, 1,1,1],  the  second  source  at  a  bearing 
of  11  degrees  was  observed  with  =  [1,  .8,  .6,  .4,  .2] .  In  this  case  the 
difference  in  intensity  may  be  viewed  as  caused  by  the  directional  pattern 
of  the  sensors  rather  than  the  directional  pattern  of  the  sources.  We 
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generated  100  independent  samples  at  a  SNR  of  20  dB.  The  spatial  spectrum, 
P(y)»  is  plotted  versus  the  angle  of  arrival  (bearing)  in  Figure  1.  Two 
very  sharp  peaks  are  observed  at  -9  degrees  and  11°  degrees.  The  associated 

estimates  of  the  intensity  vectors  are  aj  =  ( .99,1.0, .99, .99, .99)  and 

•T 

22  =  (1.0, .79, .59, .40, .20) .  The  spurious  peak  at  3  degree  is  associated 
with  8^  =  (1.0  ,0.7,0.3,— 0.1,— 0.3)  and  therefore  can  be  easily  eliminated, 
since  under  our  assumptions  the  *s  must  be  positive.  For  comparison,  we 
plotted  the  result  of  the  MUSIC  algorithm  [3]  in  Figure  2.  Since  only  one 
source  conforms  with  the  assumptions  of  MUSIC,  only  one  peak  is  observed. 
Example  2. 

Consider  Example  1  except  that  here  the  SNR  =  50  dB,  7^=11  degrees 
Y2=25  degrees  and  aj  =  22  =  (1,1, 1,1,1).  The  spatial  spectrum  is  plotted  in 
Figure  3.  We  observe  2  peaks  at  11  and  25  degrees  and  3  more  spurious 
peaks.  The  two  peaks  on  the  leftside  are  associated  with  nonphysical 
intensity  vectors  and  therefore  can  be  eliminated  by  post  processing.  The 
spurious  peak  at  18  degrees  is  associated  with  fi  =  ( .77, .94,1.0, .92, .73) 
and  therefore  is  an  ambiguous  solution.  Ambiguous  solutions  occur  whenever 
the  continuum  a(0)  ("array  manifold")  intersects  the  signal  subspace  (the 
space  generated  by  the  columns  of  A)  in  more  than  N  points  [3]. 
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Cone Ins  ions 

In  this  section  the  eigenstrnctnre  approach  has  been  used  to  obtain 
estimates  of  source  locations  as  well  as  estimates  of  the  intensity  vectors 
{Oji},  simultaneously.  The  estimates  of  may  be  useful  in  their  own 
right,  but  their  estimation  is  essential,  even  if  one  is  only  interested  in 
the  source  locations,  in  cases  where  it  is  not  appropriate  to  assume 
omnidirectionality.  For  exan^le,  whenever  a  source  is  in  the  near  field  of 
the  array,  its  radiation  pattern  can  rarely  be  assumed  omnidirectional. 
This  is  also  important  in  applications  in  which  it  is  unrealistic  to  assume 
that  the  radiation  pattern  of  each  sensor  is  accurately  known  (this  usually 
requires  frequent  calibration  and  a  large  memory). 

We  observed  that  in  some  cases  post-processing  is  required  to  eliminate 
spurious  solutions.  The  post  processing  decisions  may  rely  on  the  sign  of 
the  intensity  vectors  (8^}  and  on  any  prior  knowledge  concerning  the  source 
locations  and  expected  intensity  vectors. 
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V.  Nonuniform  Array  Processing  Via  the  Polynomial  Approach 

Recently  an  effective  technique  for  computing  the  maximum  likelihood 
(ML)  estimates  of  the  signals  was  introduced  by  Bresler  and  Macovski  [5]  and 
Eumaresan,  Scharf  and  Shaw  [17],  [18],  We  refer  to  this  technique  as  the 
"polynomial  approach*  since  it  is  based  on  expressing  the  ML  criterion  in 
terms  of  the  prediction  polynomial  of  the  noiseless  signals.  The  pol3momial 
approach  relies  on  the  assonq>tion  that  the  array  of  sensors  is  uniformly 
spaced.  It  is  well  known  [22]  that  the  optimal  sensor  configuration  is  not 
uniform  under  many  reasonable  criteria.  For  exaiiq>le,  minimum  bearing 
variance  is  obtained  by  placing  half  of  the  sensors  (with  a  spacing  of  half 
of  the  design  wavelength)  at  each  end  of  the  given  aperturej  minimum  range 
variance  is  obtained  by  placing  one  fourth-  of  the  element  at  each  end  and 
half  in  the  middle;  and  optimal  position  estimation  is  obtained  by  placing 
one  third  of  the  sensors  at  each  end  and  the  middle.  Furthermore,  when 
operating  long  uniform  arrays,  often  some  of  the  sensors  do  not  function  and 
their  outputs  must  be  ignored,  yielding  in  effect  a  sublattice  array.  In 
this  section  we  present  a  method  for  extending  the  polynomial  approach  to 
sublattice  arrays.  We  treat  the  sublattice  array  output  as  an  incon^lete 
data  observation.  Therefore  the  EM  algorithm  is  directly  applicable.  This 
algorithm  was  only  recently  applied  to  array  processing  problems  by  Feder 
and  Weinstein  [6],  However,  in  [6]  the  EM  algorithm  is  used  to  enable  the 
estimation  of  one  signal  at  a  time,  while  here  it  is  used  to  enable  the  use 
of  the  polynomial  approach  which  estimates  all  the  signals  simultaneously. 
Since  the  pol3momial  approach  is  not  widely  known,  the  basic  principles  of 
this  technique  are  briefly  reviewed  here  for  clarity.  Note  that  although  we 
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concentrate  on  the  array  problem,  all  the  results  are  equally  applicable  to 
the  corresponding  time  series  problem  discussed  in  [5],  namely,  the 
estimation  of  super in^osed  cockier  exponential  signals  in  noise. 

This  section  is  organized  as  follows.  The  pol3momial  approach  for 
processing  data  collected  over  a  uniform  array  is  described  in  V.l.  In  V.2 
it  is  shown  how  the  EM  algorithm  can  be  used  to  adapt  the  polynomial 
approach  to  the  case  of  sublattice  arrays.  Several  examples  of  our 
procedure  are  presented  in  Section  V.3,  and  Section  V.4  contains  some 
cone fusions. 

V.l  Uniform  Arrays  and  the  Polynomial  Approach 

Consider  N  narrowband  radiating  sources  observed  by  a  linear  uniform 
array  conq)osed  of  M  sensors.  The  sources  are  assumed  to  be  far  enough  from 
the  array,  con^ared  to  the  array  length  so  that  the  signal  wavefronts  are 
effectively  planar  over  the  array.  The  signal  at  the  output  of  the  m-th 
sensor  can  be  expressed  by 


N 

X  (t)  =  ^  s  {t-(m:-l)r  )  +  V  (t);  m  =  1,2, ...,M, 

m  i.  u  n  m 

n=l 

-  T/2  <  t  <  T/2, 


(31) 


where  (Sjj^(t) are  the  radiated  signals,  (■'''ni(^) additive  noise 
processes,  and  T  is  the  observation  interval.  The  delay  of  the  n-th 
wavefront  at  the  m-th  sensor,  relative  to  the  first  sensor,  is  given  by 
(m-l)Tj^.  The  parameter,  t^,  can  be  expressed  in  terms  of  the  sensor 
spacing,  d,  the  propagation  velocity,  c,  and  the  source  bearing, 
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relative  to  the  array  perpendicular  as 


X  =  (d/c)sin(y  ), 


A  convenient  separation  of  the  parameters  to  be  estimated  maybe 

obtained  by  using  Fourier  coefficients,  defined  by 


X 


m 


T 


-ja>ot 

X  (t)e  dt  . 
m 


Since  we  assume  that  the  spectrum  of  the  signals  is  concentrated  around  (Oq, 
with  a  bandwidth  that  is  small  con^ared  to  an/T,  a  single  Fourier 
coefficient  is  enough  to  con5>letely  describe  the  signals.  Taking  the 
Fourier  coefficients  of  (1)  we  obtain: 


-jw  (m-l)T 

X=>e  ^  “S+V,  m=l,2,...,M,  (32) 

m  ^  n  m 

n=l 

where  Sj^  and  Vj^  are  the  Fourier  coeffients  of  S^^Ct)  and  Vj^(t)  respectively. 
Equation  (32)  may  be  expressed  using  vector  notation  as 

X  =  AS  +  y,  (33) 


T 


where 
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5  "  ISl'  *2 . 

V.  [Vj, 

A  =  [a^^,  *2'***'-N^' 


a 

-n 


[1,  x„.  xl . x“ 

n  n  a 


n=l,2 


Nj 


X  = 
n 


In  many  cases  estimation  is  based  on  more  than  one  realization  of 
equation  (33),  corresponding  for  example  to  several  time  sanqtles  or 
observation  intervals.  In  that  case  we  use  the  index  j  to  denote  the 
different  realizations: 

Xj  =  +  Yj  j  =  1,2,...,J.  (34) 

Instead  of  estimating  directly  we  concentrate  on  estimating 

Under  the  assunq>tion  that  the  vectors  {V.jT  ^  are  i.i.d.  zero  mean  and 

"“J  J  -A 

2 

Gaussian  with  covariance  o  I,  the  maximum  likelihood  estimates  are  given  by 

J 

^^n^Ll  =  R  =  ^  lllj  -  AS.  11^  (35) 

XeUC  .  ^  J  J 

J=1 

where  ||*|j  denotes  the  Euclidean  norm  and  UC  stands  for  the  unit  circle 
which  is  the  parameter  space,  in  this  case. 
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The  minimization  required  in  (35)  is  not  trivial  since  the  vectors  {S.} 

3 

and  the  matrix  A  are  not  known  to  the  observer.  However,  whenever  A  is 
known,  R  is  minimized  by  choosing 


S.  =  (A®A)  ^A®X.  (36) 

~3  ~3 

as  the  estimate  of  Sj,  for  j=l,2,...,J,  where  (  )®  represents  the  Hermitian- 
transpose  operation.  Substituting  (36)  in  (37)  we  obtain: 

J  J 

R=  J  |jXj  -  A(A®A)"^A^J|  =  ^  (37) 

j=l  j=l 


where 


Pg  =  I  -  A(A®A)  ^A®'. 

The  polynomial  approach  relies  on  the  introduction  of  the  polynomial 
b(z)  =  bQZ^  +  b  +,..+  b|q,  whose  zeros  are  the  parameters  of  interest 
^^n^n=l*  Observe  that  by  definition  the  Mx(M-N)  Toeplitz  matrix  B  defined 


by 
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^N-1  ^0 

B®  =  '’n-1  *  •  •  ^0 

•  •  • 

•  •  « 

*^N  *^N-1  •••*  ^0 

is  orthogonal  to  A,  i.e.  B®A  =  0,  and  hence  Pg  =  B(B®B)~^B®.  Now  the 
minimization  in  (35)  can  be  expressed  in  terms  of  the  coefficients  {b£}^_Q 
as 

J 

b  =  arg  min  3  x^B(B^B^  ^B^x.,  (38) 

bee,  .  ,  ^  J 

-  b  j=l 

where  b  =  [bjjj,  bj^j, . . .  ,bo]^»  ©jj  is  the  space  of  all  the  vectors  whose 
associated  polynomials  have  zeros  only  on  the  unit  circle.  It  can  be  shown 
that  since  b(z)  has  its  roots  on  the  unit  circle,  its  coefficient  vector  is 
o-conjugate-symmetricj  i.e.  b  =  a[bo,b]^, . . .  ,bj^]®  where  a  is  a  constant  of 
unit  modulus. 

The  algorithm  for  the  minimization  required  in  (38)  is  based  on  the 
relation 

=  X  b,  (39) 

•I  «J 

where  X.  is  the  (M-N)x(N+l)  matrix  defined  by: 

*3 
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X.  =  [X.<N+1:M),  X.(N:M-1),...,X.(1;M-N)], 

J  “J  ”J  -J 

coiaponeiits 
in  (38)  we 


and  X.(k:r)  describes  a  snbvector  of  X.  consisting  of  all  of  the 

J  J 

from  the  k-th  con^onent  to  the  r-th  coniponent.  Substituting  (39) 
obtain: 


b  =  arg  min  b  C  b» 

be©, 

b 


C  = 


j=l 


(40) 


This  relation  is  used  in  the  minimization  algorithm  [5],  [16]-[18].  The 


algorithm 

starts  with  any  initial 

estimate 

b<»> 

of  b  and  proceeds  as 

follows: 

a) 

Initialization  k=0,  b=b^®^ 

b) 

Confute  C^^^  according  to 

(40)  using 

to  construct  the  matrix 

c) 

Find  =  min  b®C^^^b 

be©^ 

b 


d)  Check  convergence.  NO  -  k  =  k+1,  go  to  (b) . 

YES  -  Continued. 

e)  Find  the  roots  of  the  polynomial  b^^^^)(z)  whose  coefficients  are 
given  by  b^^'*’^^. 

In  [1]  the  relation  b  =  a[bQ,bj^, . . .  ,bjj]®  was  incorporated  in  step  (c)  to 
yield  a  simple  quadratic  minimization  problem.  We  now  turn  to  the  more 
practical  situation  of  nonuniform  arrays. 
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V.2  Snblattice  Arrays  and  the  EM  Algorithm 

In  this  paper  we  are  primarily  interested  in  the  problem  where  the 
measurements  are  taken  along  a  sublattice  arrays  of  M’  sensors.  The 
sublattice  array  may  be  described  by  a  binary  vector,  1,  of  length  M.  The 
m-th  con^onent  of  1  is  1  if  the  m-th  sensor  of  the  full  array  is  part  of  the 
subarray,  and  it  is  zero  if  the  sensor  is  missing.  Equation  (34)  may  be 
converted  to  describe  a  sublattice  array  through  a  left-multiplication  by  a 
transformation  matrix  6.  The  M'xM  matrix  6  is  constructed  by  eliminating 
all  the  zero  rows  in  diag(l).  For  example  an  array  of  three  elements  in 
positions  1,2,5  is  described  by  1  =  (1,1, 0,0,1)  and 


n 


G  = 


0 


.0 


0 

1 

0 


0 

0 

0 


0 

0 

0 


O' 

0 

1. 


Multiplying  equation  (34)  by  G  we  obtain,  for  a  given  sublattice  array,  the 
equation: 


Y.  =  GX.  =  G(AS.+V,), 
~J  "J  “J  “j 


J  —  1,2,,,,,J' 


(41) 


We  refer  to  (X.)  as  the  (unavailable)  "con^lete  data*  and  to  {Y.}  as  the 

J 

fp  rjt  T*  1' 

observed  data.  Let  Y  =  [YJ,  Y2,...,Yj]  denote  the  observation  vector  and 
0  =  [ASj^,  AS2,...,ASj]  represent  the  parameter  vector.  If  f2^(YjO)  is  the 
conditional  probability  density  function  of  y  given  6,  then  the  maximum 


likelihood  estimate  of  0  is: 
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§  =  arg  max  f  (Yjd)  =  arg  max  ln{f  (Y|§)} 


(42) 


Oae 


OsO 


where  6  is  the  parameter  space.  In  order  to  use  the  poljmomial  approach  it 
would  be  useful  to  express  Y  ia  terms  of  the  complete  data  vector 
X  =  [Xj,  relation  is  given  by 

Y  =  FX  ,  (43) 

where  F  is  the  block  diagonal  matrix  with  J  blocks: 


F 


G 


*6 


The  application  of  the  EH  algorithm  to  the  problem  at  hand  requires 

only  the  determination  of  Q(G|G').  We  recall  that  X  is  Gaussian  with  given 

2  T  T  T  T 

covariance  0*^1  and  unknown  mean  6,  where  ©  =  (©J, 

© .  =  AS .  j 
J  1 


Henc  e , 


ln{f^(X|0)}  =  -M-J’lnCno^}  --^||xj|^  +  ||e||^  -  ©^  -  X®©) . 

-  CT 


Therefore 
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Q(0|e»)  =  K  - -|(|le||^  -  -  |®e)  (44) 

a 

where  K  represents  all  the  terms  independent  of  ©  and 

I  =  EajY,©'}  =  ©'  +  F®(FF®)~^(Y-F©')  .  (45) 

Note  that  the  specific  ©  that  maximizes  (17)  is  the  same  ©  that  minimizes 
the  functional 

J  J 

“l  =  I||-»II^  =  }  lllj  -  «jll^  -  }  lllj  -  “.|l^  (46) 

3-1  j-1 

Hence,  the  M-step  of  the  EM  algorithm  maybe  performed  by  using  the 
polynomial  approach  for  minimizing  (46). 

The  following  relations  are  useful  for  the  actual  proposed  procedure. 
Using  the  block  diagonal  stimcture  of  F  and  the  relations  GG®  =  I  and  G®G  = 
diag(l),  equation  (45)  maybe  rewritten  as: 

I.  =  diag(T)©'  +  G®y.  (47) 

J  J  ^ 

where  I  is  the  con^lement  of  1  (zeros  and  ones  are  interchanged).  The 
parameter  vector  ©t  is  simply  the  estimate  of  AS-  obtained  in  the  previous 

•I  «l 

cycle  and  therefore  (47)  may  be  written  also  as: 


I!xl> 
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=  diag(l){A(A®A)"^A^.}^^^  +  =  diag(l)  {(I-B{B®B)“^B®)| .}  + 

J  J  J  J 

+  G^Yj, 

using  the  notation  of  the  poljnaoinial  approach.  As  one  would  expect, 
equation  (47)  states  that  the  coa^onents  of  X.  that  correspond  to  existing 

•I 

sensors  are  always  equal  to  the  observed  data,  i.e.  the  corresponding 
con^onents  of  Y j  . 

The  proposed  EM  algorithm  may  be  summarized  as  follows: 

(a)  Initialization:  Select  initial  values  for  find  the 

corresponding  b^®^. 

Confute:  A.=GAj  S .=(A?A2)~^A^ . j 

xj®^  =  diag(I)AS.  +  G®Y.  (See  (47)) 

1  J  J 

Set:  p=0. 

(b)  Use  the  minimization  algorithm  for  uniform  arrays: 
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(b.l)  Constuct  %.  =  [Si^^(N+l:M),...,S!^^(l:M-N)], 
J  -J  ~J 

.  ,  «  ,  (0)  ,  (p) 

set  k=0,  bj  =  • 


(b.2)  Construct  B  using  b^^ 


(k) 


J 

Coelute  C  =  }  ^Ab%-\  . 

j=l 


(b.3)  Compute  =  arg  min  b^Cb- . 

\^«b 

(b.4)  Check  convergence  of  bj^.  No  -  k=k+lj  go  to  (b.2). 

Yes  -  b^^^  =  continue. 

(c)  Construct  B  using  b^^^ 

Compute: 

|(p+l)  ^  diagdXl  -  B(BV“V)xi^^  +  G®I.. 

-J  ^  -  -J  -J 

m 

(d)  Check  the  convergence  of  Xj.  No  -  p=p+l,  go  to  (b). 

Yes  -  continue. 

(e)  Find  the  roots  of  the  polynomial  b^^^(z)  whose  coefficients  are 

given  by  b^^^ . 
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V.3  Examples 

To  illustrate  the  behavior  of  the  algorithm,  let  us  consider  two 
exan^les: 

Example  1;  Consider  a  uniform  linear  array  of  6  sensors  separated  by  half  a 
wavelength  of  the  actual  narrowband  source  signals.  Now,  assume  that  the 
two  middle  sensors  are  missing  (i.e.  =  (1  1  0  0  1  l))j  this  is  the 
optimal  configuration  for  bearing  estimation  when  the  given  aperture  is  2.5 
wavelengths  and  the  number  of  sensors  is  limited  to  4. 

The  sources  are  two  narrovband  emitters  located  in  the  far  field  of  the 
array.  One  source  is  located  at  a  bearing  of  10  degrees,  while  the  second 
source  is  located  at  a  bearing  of  25  degrees.  We  generated  only  10 
independent  samples  with  a  SNR  of  30  dB.  The  initial  guess  was  =  3°, 

=  17°.  The  algorithm  converged  to  within  one  degree  of  the  right 
result  in  8  iterations,  as  shown  in  Table  1. 

Example  2;  Consider  example  1  where  the  array  is  reconfigured  so  that  1^  = 
(101001).  Note  that  only  3  sensors  are  used  and  they  are  separated  by 
one  wavelength  and  1.5  wavelengths.  Nevertheless,  the  algorithm  converged 
to  within  one  degree  of  the  right  result  in  only  7  iterations,  as  shown  in 
Table  2.  The  initial  guess  was  =  3°,  Y2^^  ~  35°. 

V.4  Summary 

We  have  proposed  a  novel  EM  algorithm  for  the  estimation  of 
superimposed  signals  observed  by  nonuniform  arrays.  The  algorithm  is 
efficient  and  provides  accurate  results  even  when  the  number  of  samples  is 
small  and  the  sensors  are  separated  by  more  than  half  a  wavelength. 
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Note  that  convergence  theorems  exist  for  the  EM  method.  However, 
convergence  theorems  for  the  polynomial  approach  are  not  yet  available  and 
therefore  further  investigation  is  required  to  prove  the  convergence  of  the 
proposed  technique.  Finally  we  would  like  to  emphasize  that  the  EM 
algorithm  is  guaranteed  to  converge  to  a  local  maximum  of  the  likelihood 
function.  Thus  we  would  expect  that  the  algorithm  described  here  will 
converge  to  the  globally  optimum  result  only  if  the  initial  estimates  are 
good  enough.  Fast  initial  estimates  can  be  obtained  by  using  simpler 
methods  such  as  the  MOf,  MEM  or  the  MUSIC  techniques  (see  [1]  for  a  review 


of  these  methods). 


Iterations 

No. 

m 

Yi 

degrees 

^2 

degrees 

0 

3.00 

17.0 

1 

6.15 

19.38 

2 

7.29 

20.49 

3 

8.16 

21.47 

4 

8.78 

22.30 

5 

9.23 

22.95 

6 

9.55 

23.46 

7 

9.11 

23.83 

8 

9.93 

24.12 

9 

10.04 

24.32 

10 

10.12 

24.47 

T  = 


Table  1:  Evolution  of  the  algorithm  for  1 


(110  Oil) 
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Iterations 

No. 

m 

n 

degrees 

m 

Y2 

degrees 

0 

3.00 

35.00 

1 

-0.01 

18.13 

2 

3.46 

18.27 

3 

7.18 

20.16 

4 

8.74 

21.90 

5 

9.39 

23.01 

6 

9.68 

23.69 

7 

9.84 

24.10 

8 

9.92 

24.35 

9 

9.96 

24.51 

10 

9.99 

24.61 

T  ^ 


Table  2:  Evolution  of  the  algorithm  for  1 


(101001) 
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Figure  1 
Figure  2 
Figure  3 


Figure  Caption? 


Spatial  spectrum  of  the  proposed  procedure  for  two  far-field 
sources. 

Spatial  spectrum  of  the  MUSIC  procedure  for  the  case  of 
Figure  1. 

Spatial  spectrum  of  the  proposed  procedure  for  two  far-field 
sources  with  equal  intensity  vectors. 
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